#'
#' A grid based rectangle kernel smooth
#'
#' @aliases define_rectSmooth apply_rectSmooth
#'
#'@usage
#'define_rectSmooth(r,h)
#'apply_rectSmooth(rdg,dist)
#'
#'@param dist a vector of observed distances
#'@param r a vector of distance to do the kernel smooth
#'@param h length of the rectangle
#'@param rdg an rdg object which defines a specific rectangle kernel smooth
#'
#'@details
#'This is a grid-based implimentation of rectangle kernel smooth. 
#'
#'@return
#'The define_rectsmmoth returns a rdg object represents a defined rectangle smooth rule
#'The apply_rectSmooth returns rtsm contains the index of observed distances that meet the given rectangle kernel smooth at a given r position
#'
#'@examples
#' r=seq(0,98,5)
#' h=2
#' rdg=define_rectSmooth(r,h)
#' 
#' 
#' d=runif(10,0,110)
#' re=apply_rectSmooth(rdg,d)
#' 
#' plot(re,rdg,d)
#'
apply_rectSmooth=function(rdg,dist){
  #create the grid with length h
  h=attr(rdg,"h")
  r=attr(rdg,"r")
  ndg=attr(rdg,"ndg")
  
  #put the observed distances into the grid
  di=dist%/%h+1
  dist_group=lapply(1:ndg,function(x) which(di==x))
  
  #get the pottential observed distances and find the exact ones
  result=list()
  for(i in 1:length(r)){
    dindex=unlist(dist_group[rdg[i,]])
    sel=abs(dist[dindex]-r[i])<h
    result[[i]]=dindex[sel]
  }
  class(result)=c("rtsm",class(result))
  return(result)
}

define_rectSmooth <- function (r,h) {
  drgroup=seq(min(r),max(r),h)
  ndg=length(drgroup)
  nr=length(r)
  rdg=matrix(nrow=nr,ncol=3)
  for(i in 1:nr){
    rdgi=ceiling(r[i]/h)
    rdg[i,]=c(rdgi-1,rdgi,rdgi+1)
  }
  rdg[rdg>ndg]=NA
  rdg[rdg<1]=NA
  
  attr(rdg,"r")=r
  attr(rdg,"h")=h
  attr(rdg,"ndg")=ndg
  class(rdg)=c("rdg",class(rdg))
  return(rdg)
}

plot.rtsm=function(rtsm,rdg,dist){
  h=attr(rdg,"h")
  r=attr(rdg,"r")
  ndg=attr(rdg,"ndg")
  plot(x=c(1,(max(r)+2*h)),y=c(0,1),type="n",xlab="r",ylab="",main="Rectangle kernel smooth")
  segments(x0=r,y0=0,x1=r,y1=0.1)
  segments(x0=r-h,y0=0.1,x1=r+h,y1=0.1)
  segments(x0=r-h,y0=0.1,x1=r-h,y1=0.12)
  segments(x0=r+h,y0=0.1,x1=r+h,y1=0.12)
  
  points(x=dist,y=rep(0.08,length(dist)))
  num=unlist(lapply(rtsm,length))
  h=num/sum(num)
  segments(x0=r,y0=0.1,x1=r,y1=0.1+h,lwd=5)
}